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Abstract 

We present an efficient algorithm for calculating the minimum energy path (MEP) and energy 
barriers between local minima on a multidimensional potential energy surface (PES). Such paths 
play a central role in the understanding of transition pathways between metastable states. Our 
method relies on the original formulation of the string method [Phys. Rev. B 66, 052301 (2002)], i.e. 
to evolve a smooth curve along a direction normal to the curve. The algorithm works by performing 
minimization steps on hyperplanes normal to the curve. Therefore the problem of finding MEP on 
the PES is remodeled as a set of constrained minimization problems. This provides the flexibility of 
using minimization algorithms faster than the steepest descent method used in the simplified string 
method [J. Chem. Phys., 126(16), 164103 (2007)]. At the same time, it provides a more direct 
analog of the finite temperature string method. The applicability of the algorithm is demonstrated 
using various examples. 

PACS numbers: 
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I. INTRODUCTION 



The dynamics of complex systems often involve thermally activated barrier-crossing 
events that allow the system to move from one local minimum of the energy surface to 
another. At finite temperatures, the total kinetic energy accessible to the system is on the 
order of NksT, where N is the number of degrees of freedom, fee is the Boltzmann constant 
and T is the temperature of the system. However, this huge amount of energy is distributed 
over the whole system. Consequently, it fails to cross over the free energy barrier (generally 
an index-1 saddle point) and move to a different basin of attraction. A system can overcome 
a free energy barrier only when sufficient energy is localized on an activated region of the 
system. The activated region is the volume of the sample where bond breaking/formation, 
atomic re-arrangements etc. take place.- 1 ^ It is of great theoretical and practical interest to 
develop algorithms that can enable us to efficiently compute the most probable pathways 
for such transition events. For systems with relatively smooth energy landscapes, it can 
be shown that the most probable pathways are the minimum energy paths (MEP). Min- 
imum energy paths are physically relevant in the low temperature dynamics of a system 
and provides information only about the energy barrier involved in a thermally activated 
event without any consideration of the width of the channel near the saddle point or other 
entropic effects. Further, at high temperature, the energy surface becomes rugged due to 
thermal fluctuations and the presence of multiple peaks of O (fc^T) makes the concept of 
MEP irrelevant. However in such a scenario the MEP can still correspond to the path with 
the maximum likelihood. 

The problem of finding MEP and the bottlenecks for transition events can be broadly 
categorized into two classes depending on the initial conditions: (a) when only the initial 
point is known, and (b) when both the initial and final points on the energy surface are 
available. In the former case, one can resort to methods like gentlest ascent dynamics^, 
dimer method 4 , etc. to explore the energy surface. For the second category, the most 
notable examples include the string method^- and the nudged elastic band method^— . In 
this case, we are given the initial and final states of the system, and our aim is to find the 
MEP connecting these states. Since there can be multiple paths joining the end points, the 
converged MEP is dependent on the choice of the initial path. In the original string method, 



2 



a path 7 evolves as: 

^ = _Vy( 7 ) ± + rt (1) 

where, 7 is the time derivative of 7, W 1 - = W — (VV, f) f is the gradient of the potential 
perpendicular to 7, f is the unit vector parallel to the tangent to 7 and r is a Lagrange 
multiplier used to enforce a particular parametrization of the string. 

The string method is easy to implement and works well even if the initial and the final 
states are not local minimum.- This advantage has led to a class of algorithms called growing 
string method that have been tested for real systems using quantum mechanical tools.— ~— 
Further, if one is solely interested in knowing the free energy barrier and the configuration at 
the saddle point, the end of the string need not be a local minimum but some intermediate 
configuration lying in a basin other than that of the initial point. 

The string method is based on the idea of moving curves by using a steepest descent-type 
of dynamics. It should be emphasized that even though the dynamics used in the string 
method has very strong steepest decent flavor, the method does not amount to minimizing 
any energy functionals. In order to apply quasi-Newton type of ideas to accelerate the 
string method, one has to resort to the Broyden formulation for solving coupled system of 
equations rather than the BFGS type of formulation for optimization.— 

We propose a method that reduces the problem of finding MEP to an optimization prob- 
lem, so that techniques from optimization theory can now be used directly. This prescription 
is in some ways an improved version of the locally updated planes (LUP) method proposed 
earlier by Elber and co-workers.— 1 ^ In the LUP method, each intermediate configuration 
along the approximate path is relaxed by confining its motion to a hyperplane. The hyper- 
planes are selected such that they are perpendicular to the straight line joining the two end 
points. The LUP method is unstable because - (i) each intermediate configuration is relaxed 
independently and the algorithm does not prescribe any scheme to make sure that the path 
is smooth, (ii) the hyperplane selection scheme can result in the formation of kinks if the 
path has multiple local energy wells, and (iii) since there is no prescribed way to control 
the separation between the intermediate configurations, the relaxed path can have images 
clustered around the local minima. These problems severely limit the convergence to the 
correct MEP in the LUP method. 

The finite temperature string (FTS) method is like an expectation-maximization (EM) 
algorithm for curves.— 1 ^ The current version of the string method is more in line with the 
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spirit of FTS in which the sampling procedure is replaced by the minimization step. 



II. THE METHOD 

Given a curve 7 that joins two points A and B on the energy surface, let us parametrize 
7 as 7 = {X (a) : a G [0, 1]} where a is a continuous parameter. When the end points A 
and B lie in different basins of attraction, the curve 7 will pass through multiple saddle 
points that are generally the bottlenecks of transition of a system from A to B. The path 
7 is the MEP if 

(V^(X(«)) = (2) 

where (VV) ± = VV — (r, W)t is the component of VV^ normal to the path 7 and r is the 
unit tangent vector of 7. If P (a) is the hyperplane perpendicular to the path 7 at X (a), 
then (j2J) can be restated as 

X (a) = argmin V\ F(a) (3) 

This is a variational characterization of the MEP. The new version of the string method is 
based on the variational formulation in ([3]) . The initial guess path is first discretized to give 
|X°} N . At step m, the path is updated by the following procedure: 

1. Minimization: a number of minimization steps are performed on the hyperplanes 
normal to the curve at the discretization points. This gives {X!-} N . 

2. Mixing: a mixing scheme is applied to get 

X J = Xf(l-A) + AX J * (4) 
where, A G (0, 1) is the mixing co-efficient. This step helps in stabilizing the scheme. 

3. Reparametrization: redistribute the intermediate configurations along the string ac- 
cording to some metric, such as, equal spacing in configuration space, or equal spacing 
in energy space, etc. 

C _ ~\ reparametrize , 

W ► {XT 1 } (5) 

These steps are performed until convergence. Step 3 is a standard procedure in the string 
method and accounts for the displacements parallel to the curve. Steps 1 and 2, are different 
from the original and simplified string methods. 
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III. ALGORITHMIC DETAILS 



Step 1: Minimization 

The minimization procedure can be implemented in different ways. The energy of each in- 
termediate configuration can be minimized on their respective hyperplanes till convergence 
or the minimization can be performed only for few steps. Minimization algorithms with bet- 
ter convergence like, FIRE (fast inertial relaxation engine) 20 , conjugated gradient^, limited 
memory BFGS 21 , etc. can be used for this purpose. For the sake of completeness below we 
present the modified versions of BFGS and FIRE algorithms: 

BFGS: Starting with an intermediate configuration at x = X° on the path and an approx- 
imate Hessian matrix H = VW (x ) at x , the BFGS scheme involves the following steps 
until convergence is achieved :— 

(i) obtain a direction p^: H^p^ = VF(xfc), where F (x^) = — VV^(xfc). This is performed 
by obtaining the inverse of H k by applying the Sherman-Morrison scheme. 

(ii) obtain step size a k along p fc (line search) 

(Hi) update image: x fc+1 = x fc + a k p k , where p^ = p^ — (p^, Tj) Tj is the projection of the 
search direction perpendicular to the tangent (Tj) to the path at X°. 

(iv) set s k = a k p k and y fc = - [VF (x fc+ i) - VF (x fc )]. 

(v) update Hessian: 



-nfc+l = ti k + T . (bj 



yfcy fc _ H k s k slH k 

yl s k sjH k s k 

FIRE: Starting with an initial intermediate configuration at x = X° on the path, the 
following dynamical system can be used for the constrained minimization:— 



1 , * (7) 

■F-jSv + ^vF 



m 



where, Tj is the tangent to the path at X°, m is the mass and /3 is a parameter. For the 
intermediate images, the tangent Tj at X° can be approximated as 



^= IY0 +1 S"V J = 2 ^-( N - 1 ) ( 8 ) 



Step 2: Mixing 

During minimization, since each intermediate image is relaxed independently, their relax- 
ation step lengths can vary due to which the path develops kinks and is no longer smooth. 



This makes the tangent vector inaccurate leading to numerical instability. To overcome 
this difficulty, in the second step, we use a mixing scheme to have better control over step 
lengths. 

Step 3: Reparametrization 

Next, reparametrization is performed in two steps: computing the values of the parameter 
a and performing interpolation to find the reparametrized intermediate configurations. For 
parametrization by equal arc length we first obtain the length of the string: 

Lj = Lj-i + Xj - X 3 -_i (9) 

where, Lq = and j — 1,2, ...,N. The corresponding normalized parameter is then given 
by a.j = Lj/Ljy. Next, using cubic spline interpolation scheme, we obtain the intermediate 
configuration positions {X™ 4-1 } that are evenly distributed along the string according to the 
updated parameter ctj = j/N. For computational efficiency, reparametrization may not be 
enforced at each step. In fact, we have found that reparametrization is much less important 
for this version of the string method than the original version.-^ Thus, by constraining the 
minimization on the normal hyperplanes, we have already removed the strong tendency for 
the intermediate positions to move towards the local minima. So reparametrization only 
plays a minor role here for improving the accuracy by optimizing the distribution of the 
points along the curve. 

Naturally one is interested in a comparison of the performance of this version of the string 
method with the older version. The answer is that it depends on the problems one needs to 
deal with. Suffice to say that the current version provides an alternative. In addition, this 
current version of the string method is much closer in spirit to the finite temperature string 
method in which the minimization step is replaced by a sampling step.— 

IV. ILLUSTRATIVE EXAMPLES 
A. Two-dimensional potential 

To get a better understanding of the effectiveness of the selection of hyperplanes, let us 
look at the potential energy surface of a simple two-dimensional toy problem shown in Fig. 
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[2j The potential energy is given by£ 

V(x,y)={l-x 2 -y 2 ) 2 + ^- 2 (10) 
v ' x 2 + y 2 

The line AB is perpendicular to the line joining the two local minimum points (-1,0) and 
(1,0). C is at the intersection of the initial guess path and AB. If we perform a constrained 
relaxation of C along AB, there are more than one local minimum points in the vicinity. 
Consequently, during subsequent iterations the intermediate configuration can flip between 
these locally stable structures. This will make the string uneven leading to the development 
of kinks due to which it can not converge to the MEP. 

However, consider now the line EF which is normal to the path at D. Now the point 
D has only one choice to lower the energy - go from D towards E, because going towards 
F means going uphill which increases the energy. Clearly, constrained minimization along 
hyperplanes normal to the string is a better choice. After relaxation we obtain a smooth 
string joining the two local minimum points which is shown in Fig. |2] (red dashed curve). 

B. Ad-atom diffusion on (111) surface of Cu 

Next, we use the optimization-based string method to find the MEP corresponding to 
rare events taking place in nature. As an example, we study the diffusion of an ad-atom on 
the {111} surface of copper. The {111} surface is used as it has lowest surface energy in Cu. 
The inter-atomic interactions are modeled using embedded atom method (EAM) potential 
developed by Mishin et al.— 

On a pure {111} surface there are multiple sites at which the potential energy is locally 
minimum. One of them is the face center cubic (fee) hollow site and the other being the 
hexagonal close-packed (hep) hollow site. In a fee crystal, a set of three {111} planes forms 
the stacking sequence while in hep a set of two {111} planes forms the stacking sequence. 
Hence, an ad-atom occupying a surface site which corresponds to fee (hep) stacking is called 
fee (hep) hollow site. 

The simulation cell contains 512 atoms with cell dimensions of 37.57x20.45x17.71 A 3 
and cell axes parallel to [111], [110] and [112] directions. Periodic boundary conditions are 
imposed along the [110] and [112] directions and free surface conditions are imposed along 
the [111] direction. The end points of the string are configurations with ad-atoms in hep 
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hollow sites separated by about 10 A on the (111) surface. 

Fig. [3] shows the converged MEP (force norm less than 10~ 3 eV/A) as a function of the 
length of the path in the multidimensional configuration space. The MEP is obtained by 
performing constrained conjugated gradient minimization on each intermediate configura- 
tion. The red curve corresponds to the converged MEP with intermediate configurations 
evenly distributed (now shown in the figure) along the path. The blue circles denote an 
intermediate path without reparametrization. The intermediate structures in this case clus- 
ter around each other. As shown in the figure, the ad-atom in hep site is not the lowest 
energy, the structure with ad-atom occupying the fee hollow site has lower energy (~0.004 
eV). The diffusion barrier from hep hollow site to fee hollow site is 0.036 eV. Similarly, the 
diffusion barrier from fee to hep hollow site is about 0.04 eV. The results shown are without 
any mixing scheme to update the image positions, i.e. A = 1. 

C. Dislocation nucleation in a nano-wire 

Understanding the process of dislocation nucleation and multiplication is central to our 
understanding of structural stability of materials. A dislocation is a line defect which cor- 
responds to the interface between the sheared and unsheared regions of a sample. A bulk 
sample in nature generally has an astronomical number of defects, of varied dimensionality, 
present in them. In contrast, in materials of small volume the density of defects can be 
very small. Small volume materials, such as, quantum dots, nano-wires, thin films, nanos- 
tructured materials, etc. have much higher surface area (or grain boundary area) to volume 
ratio than their bulk counterparts^ 1 ^ and defect nucleation from the surface becomes an im- 
portant driving force for plastic deformation. Given the contrasting length-scales involved, 
the deformation mechanisms changes from bulk-dominated plasticity to surface-dominated 
plasticity with decrease in characteristic length scales. One such potentially important form 
of deformation in nano-wires, nano-pillars, etc. that can play a critical role in controlling 
plastic deformation is dislocation nucleation.—"— Here, we perform atomistic simulations 
of heterogeneous dislocation nucleation in a Cu nano-wire. We use a nano-wire of square 
cross-section to include the corner as a preferential nucleation site. 

The simulation cell has 28,900 atoms and the cell dimensions are 90.375x90.375x90.375 
A 3 . The cell axes are parallel to [100], [010] and [001] directions. Periodic boundary con- 
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ditions are maintained along [001] direction (i.e. wire axis) and free surface conditions are 
imposed along the remaining two directions. The interatomic interactions are modeled us- 



ing embedded atom method (EAM) potential developed by Mishin et al.— Fig. 6(a) shows 
the initial structure of the nano-wire. Central symmetry scheme of coloring is used. 28 The 
converged MEP (force norm smaller than 10 -3 eV/A) is obtained by performing constrained 
conjugated gradient minimization on each intermediate configuration along the path. 

The initial structure is a defect-free pure Cu nano-wire and the final structure is one with 
a dislocation loop (not a local minimum). The initial structure is generated by applying a 
compressive strain of 6 % along the nano-wire axis and then performing conjugated gradient 
relaxation. The final structure with a dislocation loop is generated by applying a relative 
displacement equal to a partial Burgers vector between two {111} planes. The initial guess 
path is generated using linear interpolation between the end point configurations. 

Fig. [5] shows a section of the multi- dimensional PES as a function of the distance 



between the initial point and intermediate configurations. Fig. 6(6) shows an intermediate 
configuration close to the saddle point. The activation energy barrier is 0.43 eV and the 
sheared region consists of around 30 atoms. Central symmetry scheme of coloring is used so 
only the atoms in the sheared region are visible.— More details about the thermally activated 
nature of dislocation nucleation can be found elsewhere.— 



V. CONCLUSION 

In this paper a modified version of the string method is proposed. The proposed method 
allows determining the MEP by finding the minimum energy configuration along the hyper- 
planes normal to the path. The algorithm is simple, easy to implement and provides the 
flexibility of using faster minimization algorithms. 

The applicability of the algorithm is demonstrated using a simple two-dimensional po- 
tential energy landscape. The algorithm can be easily extended to study systems of higher 
dimensions. As examples taken from nature, diffusion barriers of an ad-atom on the surface 
of Cu and the process of heterogeneous dislocation nucleation from the corner of a nano-wire 
are evaluated. For the ease of reference, codes for some sample problems are placed in a 
publicly accessed web link. 30 
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FIG. 1: Schematic representation of the optimization-based string method to find MEP. The 
intermediate configuration at X m (a) is relaxed to X m+1 on a plane P (a) which is perpendicular 
to the tangent at X m (a) . 
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FIG. 2: Contour plot of the two-dimensional potential energy landscape showing the MEP (red- 
dashed curve) and an unrelaxed path (red curve) joining the two local minimum points (-1,0) and 
(1,0). Lines AB and EF show the different possibilities of selecting hyperplanes at points C and 
D, respectively, on the string. 
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FIG. 3: Converged MEP of the ad-atom diffusion on a (111) free surface in Cu. At the end point the 
ad-atom occupies hep hollow sites that are about 10 A apart. The energy barrier for diffusion from 
hep to fee hollow site is 0.036 eV and from fee to hep hollow site is 0.04 eV. The converged MEP is 
shown in the red-curve while the blue circles show an intermediate path without parametrization. 
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FIG. 4: A snapshot of the initial structure of ad-atom on (111) surface of Cu sample. The MEP 
is shown in red dashed line. Coordination number (CN) coloring is used in this case. The atoms 
on the free surface have CN 9 while those near the ad-atom have CN 10. The ad-atom sitting on 
the hep hollow site has CN 3. 
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FIG. 5: The converged MEP of the heterogeneous dislocation nucleation from the corner of a Cu 
nano-wire. The zero temperature energy barrier is 0.43 eV for a 6 % uniaxial compressive strain 
applied along a direction parallel to the wire axis. The initial state is a defect free Cu nano-wire 
with (100) surfaces. 
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(a) (b) 

FIG. 6: Snapshots of the initial and the saddle configuration for the process of heterogeneous 
dislocation nucleation from the corner of a nano-wire. The nano-wire is under uniaxial compression 
of 6 %, the strain axis being parallel to the wire axis. The central symmetry scheme of coloring^ is 
used: 6(a) only the atoms on the surface are visible, [6(b) | the surface and the atoms in the sheared 
region, in the saddle point configuration, are visible. The viewing direction is parallel to the wire 



axis. 
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